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Abstract — We  are  developing  a  randomized  approach  to  test 
generation  for  hybrid  systems,  and  control  systems  in  general, 
inspired  by  the  Rapidly-exploring  Random  Trees  (RRTs) 
technique  from  robotic  motion  planning  which  has  proved 
successful  in  solving  high  dimensional  nonlinear  problems. 
The  approach  represents  an  automated  analysis  alternative 
for  systems  where  computing  the  reachable  set  is  intractable. 
The  standard  RRTs  method  creates  a  tree  in  the  state  space 
by  uniformly  generating  random  sampling  point  and  trying 
to  find  inputs  which  connect  them.  In  this  paper  we  propose 
a  novel  adaptive  sampling  strategy.  We  initially  bias  the 
distribution  so  that  states  near  the  “unsafe”  set  are  selected. 
We  continually  monitor  the  growth  of  the  tree.  As  the  growth 
rate  of  the  tree  declines  we  adjust  the  sampling  distribution 
to  be  less  biased.  This  adaptive  search  strategy  varies  bias 
between  “greedy”  and  global,  often  finding  test  trajectories 
more  quickly  than  the  traditional  algorithm. 

Index  Terms — Motion  Planning,  Randomized  Algorithms, 
Hybrid  Systems,  Test  Generation,  Adaptive  Biasing. 

I.  Introduction 

Hybrid  systems  provide  mathematical  models  of  em¬ 
bedded  or  software  controlled  physical  systems.  It  is  well 
known  that  the  interaction  between  discrete  and  continuous 
time  dynamics  of  such  systems  can  produce  rich  and  often 
unexpected  behavior.  For  this  reason,  as  these  systems  grow 
in  complexity  and  sophistication,  the  need  for  automated 
design  tools  increases.  The  focus  to  date  in  the  literature 
has  been  on  the  formal  verification  of  safe  operation,  via 
the  solution  of  the  reachability  problem,  initially  through 
symbolic  methods  [21],  [13]  and  later  through  numerical 
techniques  [2],  [1],  [6],  [19].  However,  the  class  of  hybrid 
systems  for  which  the  reachability  problem  is  tractable  is 
quite  limited  in  both  expressiveness  and  dimensionality. 

Test  generation  -  generation  of  the  set  of  test  inputs  to 
identify  system  faults  and  confirm  correct  system  behavior 
-  is  a  relatively  new  approach  to  analyzing  dynamic  sys¬ 
tems,  whereas  it  is  a  well  established  concept  for  software 
design  [22].  Rather  than  prove  system  safety  exhaustively, 
our  approach  is  to  try  to  generate  a  set  of  test  scenarios  that 
cause  the  system  to  satisfy  a  specification.  A  specification 
is  a  certain  property  we  intend  to  check.  For  example,  a 
specification  is  an  unsafe  set  in  safety  verification  problems. 
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We  are  developing  a  randomized  approach  to  test  genera¬ 
tion  for  hybrid  systems,  and  control  systems  in  general, 
using  techniques  from  robotic  motion  planning  which  have 
proved  successful  in  solving  high  dimensional  nonlinear 
problems.  The  merit  of  this  approach  as  compared  with 
exhaustively  proving  safety  through  reachability  analysis 
is  that  decidability  issues  do  not  come  into  play  because 
we  are  not  attempting  to  represent  the  entire  reachable  set. 
The  drawback  of  the  approach  is  that  it  is  a  semi-decision 
method,  meaning  that  we  can  only  disprove  system  safety 
by  counter  example  -  safety  cannot  be  proved.  Despite 
this  drawback  randomized  approaches  hold  great  promise 
for  addressing  complex  nonlinear  real-world  problems  for 
which  trial  and  error  testing  is  not  sufficient  and  formal 
analysis  is  intractable. 

The  basic  algorithm  introduced  here  is  inspired  by  the 
Rapidly-exploring  Random  Trees  (RRTs)  [15],  [16]  tech¬ 
nique  from  robotic  motion  planning.  The  standard  RRT 
method  creates  a  tree  in  the  state  space  by  uniformly 
generating  random  sampling  points  and  trying  to  find  inputs 
which  connect  them.  Other  sampling  strategies  which  bias 
the  samples  in  a  region  closer  to  the  goal  state  have  been 
tried  in  [17]  and  [5]  with  some  success.  However  it  is 
difficult  to  decide  on  the  biasing  a  priori.  In  many  problems 
a  “greedy”  goal  biased  approach  will  rapidly  yield  the 
most  obvious  solution.  However,  in  problems  with  dynamics 
which  are  not  small  time  locally  controllable  or  state  spaces 
which  are  not  simply  connected  a  more  global  sampling 
strategy  is  required  to  discover  a  solution  trajectory. 

The  primary  contribution  of  this  paper  is  to  propose  a 
new  adaptive  sampling  algorithm  (Section  IV  and  V)  which 
adjusts  the  sampling  distribution  between  heavily  biased 
toward  the  “unsafe”  region  (i.e.,  greedy)  and  the  traditional 
uniform  distribution  (i.e.,  global)  depending  on  the  growth 
rate  of  the  tree.  The  outline  is  as  follows.  In  Section  II, 
we  formally  define  the  test  generation  problem  for  hybrid 
systems.  In  Section  III  we  point  out  the  connection  to  robot 
motion  planning  and  review  the  RRT  algorithm  as  well  as 
our  previous  work  in  the  area,  especially  with  regard  to 
estimating  coverage.  Section  IV  and  V  describe  the  new 
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adaptive  sampling  algorithm-  the  primary  contribution  of 
the  paper.  In  Section  VI  some  computational  statistics  for 
the  algorithm  are  presented  for  two  example  problems:  the 
thermostat  problem  from  the  hybrid  system  literature  and  a 
problem  involving  an  unmanned  aerial  vehicle. 

II.  Problem  Statement 

Definition  2.1:  We  define  a  Finite  Time  Hybrid  Con¬ 
trol  System  (modified  from  the  Hybrid  Automata,  see  [18]) 
as  a  tuple  H  =  {X,  Q,  U,  T,  Init,  /,  Inv,  E,  G, R )  where 

•  X  C  is  a  set  of  continuous  variables; 

•  Q  C  N  is  a  set  of  discrete  variables  which  index  the 
system  modes; 

•  U  C  K171  is  a  compact  set  of  possible  input  values; 

•  T  =  [to,  tf]  C  R  is  a  compact  time  interval  the  system 
evolves  over; 

•  Init  C  X  x  Q  is  a  set  of  possible  initial  conditions; 

•  /  :  X  xQ  xU  — >  M.N  is  a  vector  field  which  prescribes 
the  time  derivative  of  the  continuous  variables  (i.e., 

x  =  f{x,q,u))\ 

•  Inv  :  Q  — >  2X  assigns  to  each  q  £  Q  an  invariant  set; 

•  E  C  Q  x  Q  is  a  collection  of  edges  describing  the 
possible  discrete  transition  (a.k.a.-  mode  switches); 

•  G  :  E  — >  2X  assigns  to  each  e  =  (q,  q')  £  E  a  guard; 
and 

•  R  :  E  x  X  — *  2a  assigns  to  each  e  =  ( q ,  q')  £  E  a 

reset  relation. 

Throughout  this  paper  we  refer  to  (x,  q)  as  the  state  of  the 
hybrid  system.  Note  that  we  use  the  term  “input  signal”  in 
the  most  general  sense  in  that  it  can  include  yet  unspecified 
feedback  control  inputs,  human  in  the  loop  type  inputs, 
disturbances,  etc. 

Problem  2.2:  Testing  Problem:  Given  a  tuple 

(H,x°,q°,S),  where 

•  H  =  (X,Q,U,T,  Init,  f,  Inv,  E,G,  R)  is  a  finite 
time  hybrid  control  system, 

•  a;0,  q°  £  Init ,  and 

•  /S'  is  a  specification  set, 

the  goal  is  to  determine  an  open  loop  control  law  U  :  T  — > 
U  such  that  3t  £  T  for  which  ( x(t),q )  £  S. 

In  other  words,  the  goal  is  to  determine  a  counter-example 
-  an  input  sequence  which  will  cause  the  system  to  fail 
by  entering  S  -  if  one  exists.  However,  in  order  to  make 
the  problem  algorithmically  tractable,  instead  of  searching 
the  set  of  all  possible  functions  U  :  T  — ►  U,  the  search 
must  be  restricted  to  some  subset  of  functions  with  finite 
dimensional  parametrization. 

For  the  sake  of  convenience  we  make  three  additional 
assumptions.  First,  assume  X  x  Q  C  R^  x  N  is  defined  in 
such  a  way  that  a  point  in  M.N  x  N  can  be  easily  tested  for 
membership  in  X  xQ.  Second,  assume  the  specification  set 
S  can  be  defined  as  the  sub-level  set  of  some  function  S  = 
{(a:,  q)\x  £  X ,  q  £  Q ,  s(x,  q)  <  0}.  Finally,  we  restrict  our 
search  over  U  to  piecewise  constant  functions  of  time  with 
k  segments,  each  of  time  duration  Si.  Thus,  instead  of  the 


continuous  map  U ,  we  consider  the  search  over  U  :  T  — >  U, 
as  the  search  for  a  k-vector  of  parameters.  With  ul  £  U 

u  =  [u\u2, . . .  ,uk]T 
so  the  input  u(t)  is  given  by 

u(t)  =  ul  £  U  if  to  +  (i  —  1  )St  <t<t.0  +  ( i)St 
for  *  =  1, . . . ,  k. 

III.  Background  and  related  work 
A.  Motion  planning  and  the  RRT 

Motivated  by  the  similarities  between  the  Testing  Prob¬ 
lem  for  hybrid  systems  introduced  in  Section  II  and  the 
motion  planning  problem  from  the  robotics  literature  (in 
particular  see  [7]),  we  apply  the  Rapidly-exploring  Random 
Trees  (RRTs)  algorithm  [15],  [16]  to  the  testing  problem. 
Figure  1  illustrates  the  concepts  and  a  very  basic  algorithm 
is  given  in  Algorithms  1  and  2.  Note  that  p  is  a  suitable 
metric  function;  and  the  notation  ( x,q)  +  f  f(x,q,u)dt 
means:  using  ( x,q )  as  an  initial  condition,  simulate  the 
evolution  of  the  system  for  6t  seconds  using  u(t)  £  U  as 
the  control  input.  Various  versions  of  the  algorithm  can  be 
generated  using  different  metrics  ( p )  or  random  distributions 
(referred  to  below  as  pdf). 

The  benefits  of  the  RRT  algorithm  are  as  follows. 

•  Because  the  algorithm  works  directly  in  the  space  of 
admissible  inputs,  it  is  applicable  to  any  type  of  control 
system,  unlike  other  algorithms  [21]  which  are  only 
useful  for  systems  with  a  specific  structure. 

•  It  has  been  shown  [15]  to  be  probabilistically  complete, 
meaning  that,  if  the  problem  is  solvable,  the  probability 
of  the  algorithm  finding  a  solution  approaches  one  as 
the  number  of  nodes  approaches  infinity. 

•  Unlike  testing  schemes  which  use  pure  randomization 
of  the  inputs  the  RRT  algorithm  exhibits  a  Voronoi 
Bias  [16]  -  meaning  that  it  quickly  visits  unexplored 
regions  of  the  state  space. 

This  algorithm  has  experienced  widespread  success  in  solv¬ 
ing  a  variety  of  high  dimensional  and  nonlinear  problems  in 
motion  planning  and  has  recently  been  applied  to  controller 
synthesis  problems  for  hybrid  systems. 


Algorithm  1  Grow  Test  Set  T 

Initialize  RRT:  T.addVertex(a;0,  q°) 
while  fl(x,q)  £  T  such  that  s(x,q)  <  0  do 

ExtendfT,  pdf) 

end  while 


B.  Coverage  measures 

Because  many  testing  problems  have  no  solution  we 
must  decide  when  to  terminate  Algorithm  1.  One  way 
to  decide  when  to  terminate  the  algorithm  is  based  on 
how  well  the  tree  covers  the  state  space.  We  review  the 
concept  of  coverage  and  developments  from  [8]  here.  The 
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Fig.  1.  The  Testing  algorithm  (inspired  by  the  RRT  [15],  [16]).  The  test 
set  is  represented  as  a  tree  T  with  nodes  as  states  ( x ,  q )  and  edges  as 
inputs  u  G  U.  First  a  new  state  is  generated  at  random,  (xrand,  qrand) 
( top  left )  from  a  given  probability  density  function  (denote  as  (pdf)). 
The  algorithm  then  determines  the  closest  state,  (xnear ,  qnear)  in  the 
tree  to  the  random  state  ( top  right).  It  determines  which  u  E  U  brings 
(Xnear,qnear)  closest  to  (xrandqrand)  (bottom  [gfty  nnew  is  applj^ 
for  a  duration  5t  and  the  new  node  (xnew ,  qnew)  and  edge  nnew  are 
added  to  the  tree  ( bottom  right). 


Algorithm  2  Extend(T,  pdf) 


™ rand  nrand 
x  ,  q 

™ near  nnear 
x  ,  q 

nnew  _  argminue(7 


/  f(x,q,u)dt)} 

( xnew,qnew )  =  (x 

T.addVertex(xnew: 
T.addEdge(nnew,  (x 


random(pdf) 

nearestNeighbor(T,  (x' q' 

{p{  (; xrand  qrand ^ 


*)) 


gnear j  | 


near  ^near 


)  +  JSt  f(x,q,nnew(t))dt 
)  -c  ( xnew,qnew )) 


near  nnear 
•>  H 


Discrepancy  (a  concept  from  the  Monte  Carlo  literature)  is 
mentioned  in  [14]  but  it  is  too  expensive  to  compute  online. 
Another  appealing  idea  to  measure  the  growth  or  coverage 
of  the  RRTs  is  to  compute  the  volume  of  the  convex  hull. 
Unfortunately  the  convex  hull  is  more  indicative  of  the 
disttibution  of  the  points  than  it  is  of  the  coverage.  In  [17]  a 
variant  of  the  convex  hull  is  explored.  Rather  than  compute 
the  hull  of  all  tree  vertices,  vertices  are  grouped  according 
to  their  depth  from  parent  nodes.  The  union  of  these 
hulls  clearly  provides  a  better  approximation.  However  the 
selection  of  the  grouping  is  somewhat  arbitrary.  It  is  not 
clear  how  to  relate  the  union  to  coverage  due  to  possible 
overlaps.  Dispersion  measures  the  largest  unexplored  region 
(see  [10]  or  more  recently  [14]),  which  was  introduced  in 
the  Monte  Carlo  literature.  Assuming  we  have  a  sample  set 
X,  which  contains  N  points,  over  the  space  X,  it  is  defined 
as 

p(X,  p)  =  sup  min  p{x,  x)  (III.  1) 

xex  sex 

and  can  be  thought  of  as  the  radius  of  the  largest  empty 
ball  whose  center  lies  in  X.  We  reject  it  for  computation 
on  two  grounds:  (1)  it  is  impractical  to  compute  in  high 
dimensions;  and,  (2)  it  is  an  overly  conservative  coverage 
measure  because  it  only  considers  the  largest  ball.  For 
example,  in  Figure  2  the  left  and  right  panels  represent  two 
sample  sets  with  the  same  dispersion.  Obviously  the  sample 
shown  on  the  right  covers  the  state  better. 

In  [8]  we  introduced  a  new  coverage  measure  which  is 
used  in  this  paper  as  termination  criteria  for  the  Testing 


Fig.  2.  Two  sample  sets  which  have  the  same  dispersion  (the  radius  of 
the  largest  empty  ball  drawn  with  dashed  line).  The  set  on  the  left  clearly 
has  inferior  coverage  to  the  set  on  right. 


Fig.  3.  A  grid  is  superimposed  on  the  state  space.  The  shaded  regions 
indicate  unreachable  sets.  The  distances  from  the  grid  points  to  the  closest 
nodes  are  dj  (shown  as  dashed  arrows)  and  the  grid  spacing  is  S. 


Algorithm.  We  begin  by  overlaying  a  grid  containing  ng 
points  with  spacing  6  on  the  state  space  (see  Fig.  3).  We 
calculate  the  minimum  distance  from  each  grid  point  j  to 
the  set  of  nodes  in  the  tree,  dj.  The  quantity  mirifciy, S) 
may  be  thought  of  as  the  radius  of  the  largest  ball  centered 
at  each  grid  point  which  does  not  contain  a  ttee  node  or 
another  grid  point.  Given  a  ttee  T  we  define  its  coverage 
at  the  kth  iteration  of  Algorithm  1  as 

d  i  n9 

which  is  the  average  of  all  the  node  distances,  normalized  by 
the  grid  spacing.  Our  measure  is  similar  to  an  approximation 
of  an  “average”  dispersion  (see  equation  (III.  1)),  but  far  less 
conservative  and  faster  to  compute.  We  can  set  a  threshold, 
such  that  if  c(7fe)  <  c  the  algorithm  should  terminate.  If 
there  are  unreachable  regions  larger  than  a  grid  cell,  c(T) 
will  approach  a  non-zero  constant  as  k  — >  oo.  In  such  a 
case  we  look  at  the  change  in  c  over  the  previous  nc  nodes 

A cfc  =  —  (c(7fc)  -  c(Tk-nc))/nc.  (HI. 3) 

And  define  a  threshold  such  that  if  A cb  <  A  the  algorithm 
should  terminate.  An  exact  relationships  between  c  and  the 
true  coverage  is  unknown  at  present;  but  it  would  likely 
require  computing  the  reachable  set  which  is,  of  course, 
assumed  to  be  unknown.  Overall  one  of  the  advantages 
of  this  measure  is  that,  the  grid  size  can  be  as  fine  or 
coarse  as  one  chooses.  Finer  grids  will  tend  to  require 
more  distance  queries  but  are  more  accurate  indications  of 
coverage.  Of  course  grids  can  be  generated  in  the  “output” 
or  specification  space  to  measure  coverage  there  as  well. 


1168 


C.  Related  work 


Aside  from  general  work  on  hybrid  system  verification 
and  software  testing  mentioned  earlier,  the  idea  of  using  test 
generation,  especially  RRTs,  to  analyze  hybrid  systems  is 
new.  In  [20]  a  genetic  algorithm  approach  is  used.  The  first 
published  work  using  RRTs  for  analyzing  hybrid  systems 
is  [5],  [17].  In  a  similar  vein,  a  blimp  system  control 
law  was  validated  under  unpredictable  but  bounded  distur¬ 
bances  [12],  In  [4],  the  reachable  set  for  aircraft  collision 
avoidance  problem  was  obtained  and  several  extensions  of 
the  RRT  approach  were  mentioned.  [3]  applies  the  method 
to  biological  networks  and  [8]  is  the  first  work  to  consider 
synthesizing  time  invariant  parameters  in  addition  to  inputs. 
It  also  is  the  first  work  to  address  the  issue  of  estimating 
coverage  and  termination  criteria.  However  all  of  these 
works  use  a  fixed  sample  distribution. 

IV.  Sample  Biasing 

Traditionally  randomized  motion  planning  algorithms  use 
a  uniform  sample  distribution  to  generate  xrand  which 
tends  to  grow  the  tree  in  such  a  way  that  the  entire  state 
space  is  covered.  In  robot  motion  planning  applications,  the 
presence  of  obstacles  often  precludes  “greedy”  solutions 
which  involve  the  robot  proceeding  directly  to  the  goal. 
The  use  of  the  uniform  distribution  increases  the  algorithm’s 
robustness  when  solving  problem  in  which  the  state  space 
is  not  convex  or  the  system’s  dynamics  are  not  small  time 
locally  controllable  (STLC).  However,  when  the  state  space 
is  convex  and  the  system  is  small  time  locally  control¬ 
lable,  using  a  uniform  distribution  entails  much  wasted 
computation  because  the  “obvious”  solution  is  not  explored 
immediately. 

Since  the  goal  in  testing  is  to  find  trajectories  which  enter 
the  unsafe  region  and,  unlike  the  motion  planning  problem, 
we  assume  the  state  space  is  simply  connected,  another 
option  would  be  to  use  a  distribution  which  was  biased 
in  such  a  way  to  generate  the  majority  of  its  samples  in 
the  regions  defined  by  s(x,q)  <  0  (in  robotics  this  would 
be  akin  to  the  goal  region).  The  possible  drawback  of  a 
biased  search  scheme  (or  any  greedy  search  strategy)  is  that 
it  can  get  stuck  in  “local  minima”  for  non-small  time  locally 
controllable  systems.  Since  the  notion  of  small  time  local 
controllability  may  be  state,  mode,  or  system  dependent 
and  is  difficult  to  assess  the  utility  of  biasing  a  priori,  we 
elect  to  use  a  parameterized  distribution  whose  bias  can  be 
easily  altered.  Our  probability  density  function  B(x;p,/3), 
resembles  a  Gaussian  over  a  compact  domain,  a  <  x  < 
b.  Define  IV(x;  p,  cr(/3))  as  the  Gaussian  distribution  with 
mean  p  and  standard  deviation  o.  Define  the  probability 
density  function 

(  N(x;p,<t{/3))+ 

B{x-,p,j3)=l  Ct/(b  —  a),  a<x<b  (IV.  1) 

|  0  else 


Fig.  4.  The  distribution  B(x\  ji,  f3)  with  //  =  0  and  various  values  of  a. 
Note  that  the  function  has  compact  support. 

The  last  term,  Ct/ib  —  a),  is  added  to  ensure  that  the  area 
under  the  curve  is  equal  to  one.  Ct  represents  the  area  of 
the  truncated  portions  above  x  =  b  and  below  x  =  a. 

/a  /*oo 

N(x;  p,a)dx  +  /  N(x)  p,a)dx. 

-oo  J  b 

Obviously  p  should  be  selected  so  that  s(p,q)  <  0. 
In  Section  V  we  discuss  how  to  select  o.  Note  that 
by  changing  0  <  o  <  oo  we  can  effectively  alter  the 
distribution  to  be  completely  deterministic  ( xrand  =  p), 
completely  uniform,  or  somewhere  in  between.  Figure  4 
illustrates  the  shape  of  B(x;p,/3)  with  different  values  of 
o.  Random  states  can  be  generated  according  to  equation 
(IV.  1)  by  rejecting  random  points  generated  outside  the 
compact  domain  in  conjunction  with  any  standard  random 
normal  generator  in  each  dimension.  Note  that  qrand  can 
be  generated  using  a  similar  distribution  and  rounding  the 
result  to  the  nearest  integer. 

V.  Adaptive  Sampling  Algorithm 

An  adaptive  sampling  algorithm  should  begin  with  a 
biased  search  scheme  and  then,  if  the  tree  is  growing 
rapidly,  maintain  or  increase  the  bias;  if  the  growth  rate 
slows,  a  more  uniform  sampling  strategy  should  be  used. 
To  adjust  the  bias  we  alter  o  in  equation  (IV.  1)  using 

&  (1  /-^)(^max  ^min)  T  ^mim  (V.  1) 

where  ermax  and  <rmin  are  user-defined  values  of  the  max¬ 
imum  and  minimum  standard  deviations  and  (3  £  [0, 1]  is 
the  biasing  factor  which  we  define  below. 

Intuitively,  a  large  bias  is  effective  in  the  RRT  algorithm 
when  it  is  possible  to  steer  the  system  to  the  set  such  that 
s(x,q)  <  0  while  low  bias  is  more  advantageous  when  it 
is  difficult  to  do  so.  So  we  must  gauge  the  ease  with  which 
the  system  can  be  steered  to  a  given  state.  We  measure 
how  well  biasing  works  by  monitoring  the  growth  of  tree 
when  the  random  state  is  in  the  set  satisfying  s(x,q)  <  0. 
Consider  Figure  5.  The  vector  (xrand  —  xnear)  represents 
the  desired  direction  of  tree  growth  in  a  given  iteration  of 
Algorithm  1,  while  (xnew  —xnear)  is  the  actual  direction  of 
growth.  Therefore  if  the  absolute  value  of  the  angle  between 
these  two  vectors,  \0\  is  small  the  system  is  easily  steered  to 
the  random  state;  while  \9\  >  7r/2  means  that  the  tree  is  not 
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growing  in  the  desired  manner.  This  quantity  is  convenient 
because  it  is  naturally  bounded,  [0,7r],  and  unitless.  We 
update  the  amount  of  biasing  for  every  Ns  iterations  of 
the  RRT  algorithm,  where  Ns  is  user  defined  number.  For 
each  of  Ns  iterations,  we  can  compute  the  average  of  the 
absolute  value  of  6  over  the  n0  iterations  where  random 
states  are  generated  inside  the  set  defined  by  s(x,q)  <  0 
and  use  this  to  compute  the  bias 

P  =  tt/2 -min((l/^)X]j!fi  N^/2)  (y2) 

7r/2 

where  9i,...,6nf}  are  angles  measured  at  each  iteration. 
However,  for  systems  whose  state  does  not  evolve  over  the 
Euclidean  space,  there  is  no  meaningful  way  to  compute 
such  an  angle.  A  more  general  notion  of  measuring  the 
ease  of  steering  the  system  to  a  given  state  is  depicted 
in  Figure  6.  If  in  a  given  iteration  p(xnear ,  xrand)  > 
p{xnew ,  xrand),  where  p  is  a  metric  function,  we  call  such 
an  iteration  successful  because  the  tree  has  grown  toward 
xrand.  We  count  the  number  of  successful  iterations  ns,  out 
of  the  rip  iterations  and  compute 

0=—  (V.3) 

n0 

If  xnear  is  unsuccessful  or  the  best  candidate  for  xnew  from 
xnear  js  ajreacjy  in  a  tree  in  the  above  test,  we  eliminate  the 
xnear  from  consideration  as  xnear  for  the  testing  in  future 
iterations  to  prevent  it  from  being  chosen  repeatedly. 

The  appropriate  definition  of  (3  can  be  used  in  Algo¬ 
rithm  3  to  adaptively  alter  the  bias. 


Algorithm  3  Grow  Test  Set  T 
Initialize  RRT:  T.addVertex(a;0,  q°) 

c  =  oo,  Ac  =  oo,  0  =  1 

while  €  T  s.t.  s(x,q)  <0)Ac>cAAc>A 

do 

(J  —  (1  /?)(<Tmax  crmjn)  T  Cmin, 

pdf  =  B(x ;  p,  /3) 

ExtendfT,  pdf) 

Compute  c,  Ac 
if  new  iteration  set  then 
Update  /3 

end  if 
end  while 

Return  T 


VI.  Examples 

In  this  section  we  solve  two  example  problems  and  exam¬ 
ine  the  computational  efficiency  of  the  adaptive  sampling  al¬ 
gorithm  compared  to  fixed  biasing  approaches.  Specifically 
each  example  problem  is  solved  10  times  using  Algorithm  3. 
For  comparison  purposes  the  same  problems  are  also  solved 
10  times  using  a  variety  of  traditional,  fixed  sampling 
distribution  RRT  methods  such  as  the  one  in  Algorithm  1. 
The  fixed  bias  solutions  are  computed  using  the  distribution 


X“” 

Fig.  5.  If  |0|  is  small  on  average  (close  to  zero)  it  is  an  indication  that 
the  system  can  be  steered  in  the  desired  direction  and  a  biasing  strategy 
is  appropriate.  If  \0\  >  ir/2  it  is  an  indication  that  the  system  cannot  be 
steered  in  the  desired  direction  and  a  more  uniform  sampling  strategy  is 
appropriate. 


Fig.  6.  A  generalization  of  Figure  5  involves  comparing  two  distances 
to  determine  if  the  tree  is  successfully  growing.  If  p(xrlear,  xrand)  > 
p(xnew ,  xrand)  the  iteration  is  considered  successful. 

in  equation  (IV.  1),  with  er  =  {1(6  —  a),  3(6  —  a),  6(6  — a)}  - 
heavy  bias,  medium  bias,  and  a  nearly  uniform  distribution. 
The  fixed  bias  approach  is  similar  to  [17].  The  average 
number  of  tree  nodes  and  total  computation  time  on  a 
1GHz  PC  are  tabulated.  In  each  problem,  we  set  p  so  that 
it  was  in  the  center  of  the  region  S,  defined  by  s(x,  q)  <  0. 
U  was  discretized  into  100  possible  values  and  the  nnew 
was  computed  by  exhaustive  search.  The  algorithm  in  [9] 
was  used  to  detect  states  in  S.  All  trials  converged  to  a 
solution. 

A.  Example  1:  the  thermostat 

The  hybrid  automata  model  of  a  thermostat  has  been  a 
popular  example  in  the  verification  literature  [11],  Figure  7 
shows  the  system  model.  The  discrete  state  space  consists  of 
two  modes  Q  =  {1,2}  which  denote  the  heater  being  “on” 
or  “off  respectively.  Icl3  where  X\  is  the  temperature 
in  the  room,  X2  is  the  elapsed  time,  and  a; 3  is  a  timer  that 
measures  the  cumulative  amount  of  time  the  heater  has  been 
on  for.  U  consists  of  uon  =  [2,4];  and  u0ff  =  [-3,-1], 
The  values  uon  and  u0f  /  represent  the  possible  heating  and 
cooling  rates  in  the  two  modes.  They  are  unknown  due 
to  unmodeled  disturbances.  The  guards  gei  :  X\  <  1  and 
<?e2  :  X\  >  3  are  associated  with  the  edges  e\  :  2  — >  1  and 
e2  :  1  — >  2  respectively.  They  represent  the  temperature 
at  which  the  heater  should  turn  on  or  turn  off.  In  [11]  a 
symbolic  verification  tool  is  used  to  answer  the  question: 
“After  an  initialization  period  of  two  minutes,  is  it  possible 
for  the  heater  to  be  on  for  more  than  two  thirds  of  the  total 
time  at  any  point  dining  the  first  hour  of  operation?”  Such 
a  question  may  be  important  from  an  energy  consumption 
point  of  view.  Therefore  the  specification  set  S,  defined  by 
s(x,  q )  <  0,  is 

S  =  {si  =  x  £  X\  —  2/3x2  +  X3  >  0  A  S2  =  X2  —  2  >  0}. 

Aside  from  being  a  classical  verification  example,  the 
scenario  is  interesting  for  two  reasons.  First,  the  system 
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Fig.  7.  The  thermostat  hybrid  automata. 


Fig.  8.  The  solution  of  the  thermostat  counter  example  via  the  randomized 
algorithm. 

has  quite  nontrivial  dynamics,  since  the  discrete  dynamics 
play  a  significant  role  in  system  behavior.  For  example  the 
choice  of  uon  or  u0ff  only  influences  Xi  which  does  not 
explicitly  appear  in  S.  It  is  only  by  influencing  the  switching 
behavior  that  the  inputs  change  S.  Second,  the  set  of  all 
“bad”  trajectories  only  intersects  the  failure  region  by  a  very 
small  margin,  making  the  hunt  for  a  counter  example  quite 
difficult. 

Note  that  the  true  solution  lies  in  the  function  space  of 
piecewise  constant  inputs,  so  our  approximation  of  U  in  no 
way  impacts  the  solution.  The  test  generator  algorithm  suc¬ 
cessfully  computed  a  counter  example  as  seen  in  Figure  8. 
The  initial  conditions  were  q°  =  1  (on),  x°  =  [2  0  0]T. 
We  use  Ns  =  30  and  nc  =  30  to  ensure  that  the  means 
and  differences  were  statistically  significant,  with  c  =  0.01 
and  A  =  0.01.  The  Euclidian  distance  was  used  for  metric. 
/ 3  is  computed  via  (V.2).  Table  I  shows  the  computational 
statistics  for  the  various  distributions.  The  adaptive  strategy 
is  slightly  better  than  the  fixed  heavy  bias  and  far  superior  to 
the  other  fixed  bias  distributions  by  a  factor  of  two  or  more. 
It  is  important  to  stress  that  there  is  no  way  to  know  which 
bias  level  is  most  effective  a  priori  for  a  given  problem  so 
the  adaptive  method  is  most  appropriate. 

B.  Example  2:  hovercraft 

For  our  second  example  we  consider  a  scenario  in  which 
we  would  like  to  determine  if  a  hovercraft  (see  Figure 
9)  under  high  wind  conditions  reaches  some  “bad”  region 
such  as  a  turbulent  area.  We  denote  it  by  goal  region.  We 
assume  constant  altitude  flight  so  the  craft  has  6  states, 
x  =  (xi,X2,0,Vi,V2,u>)  and  the  dynamic  equations  are 

mi)  1  =  (/i  +  f2)  COS (9)  +  fXlair{x,  Vair{x)) 

rni)  2  =  (/l  +/2)sin(0)  +  fxzGbiript)  ^air  (*)) 

Ju  =  {f 2  -  fl)l  +  Tair(x,  Vair(x)) 


Sampling  Method 

No.  of  Nodes 

Computation  time  (sec) 

Uniform 

356 

126.4 

Med.  Bias 

164 

67.3 

Heavy  Bias 

112 

40.4 

Adaptive  Bias 

88 

37.1 

TABLE  I 


Thermostat  Example:  A  comparison  of  the  sampling 

STRATEGY  INTRODUCED  HERE  (ADAPTIVE  BIAS)  TO  OTHER 
FIXED-BIAS  SAMPLING  STRATEGIES,  AVERAGED  OVER  10  TRIALS  ON  A 
1GHz  PC. 


The  control  inputs  are  u  =  [fi  f2]T  (forward  actuating 
forces)  and  U  =  [—10, 10]  x  [—10, 10].  Forces  due  to  wind 
disturbances  in  the  x\,  x2  and  9  directions  are  fXiair , 
fX2air ,  and  Tair  whose  exact  expressions  are  omitted  for 
brevity  but  are  quadratic  in  the  difference  between  the 
craft’s  velocity  and  the  wind  velocity  valr  and  vary  with 
the  orientation  of  the  craft.  Note  that  the  state  is  partitioned 
into  two  regions  (indoor  and  outdoor)  which  define  the  wind 
velocity  differently: 

_  f  [—avx 2  (3vxi]T,  \/(xi)2  +  (x2)2  <  100 
Vair  [0  0]r,  ^(xi)2  +  (x2)2  >  100  ' 

The  biasing  factor  (3  is  computed  via  (V.3).  The  initial 
state  is  [0  0  0  0  0  0]T  and  the  goal  region,  which  is  the 
specification  set,  is  x\  €  [190,  200],  x2  €  [0,  10].  Ns  =  50 
and  nc  =  50  are  used  with  c  =  0.01  and  A  =  0.01. 

Figure  10  shows  the  solutions  of  the  problem  with 
the  uniform  sampling  distribution  and  adaptive  algorithm. 
Figure  11  shows  how  (3  changes  as  the  algorithm  evolves. 
Intuitively  biasing  is  not  effective  in  mode  1  when  the  craft 
is  subject  to  high  wind  forces  but  it  is  very  effective  when 
there  is  no  wind  in  mode  2.  The  adaptive  algorithm  is  able 
to  exploit  the  situations  in  which  biasing  is  effective.  As 
shown  in  Table  II,  adaptive  biasing  algorithm  improves  the 
efficiency  of  RRT  method  compared  to  other  bias  strategies 
rather  dramatically. 


Sampling  Method 

No.  of  Nodes 

Computation  time  (sec) 

Uniform 

3556 

1753.5 

Med.  Bias 

1017 

490.2 

Heavy  Bias 

912 

408.3 

Adaptive  Bias 

678 

342.5 

TABLE  II 


Hovercraft  Example:  A  comparison  of  the  sampling 

STRATEGY  INTRODUCED  HERE  (ADAPTIVE  BIAS)  TO  FIXED-BIAS 
SAMPLING  STRATEGIES,  AVERAGED  OVER  10  TRIALS  ON  A  1GHZ  PC. 
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Fig.  10.  RRTs  of  the  hovercraft  problem  with  uniform  sampling  (left) 
and  with  adaptive  bias  (right). 


Fig.  1 1 .  The  evolution  of  the  bias  (3  for  the  hovercraft  problem. 


VII.  CONCLUSION 

In  this  paper  we  introduce  the  idea  of  test  generation 
for  hybrid  systems,  and  control  systems  in  general,  as  a 
practical  alternative  for  scenarios  in  which  verification  via 
reachability  is  not  practical.  The  goal  is  to  search  for  a 
single  unsafe  trajectory  -  a  counter  example  to  the  safety  or 
performance  specification.  The  Rapidly-exploring  Random 
Trees  (RRTs),  a  well  established  algorithm  in  robot  motion 
planning,  was  applied  to  a  classical  verification  problem 
with  considerable  success.  Two  novel  refinements  to  the 
traditional  incarnation  of  the  RRT  algorithm  are  introduced 
to  adapt  the  method  to  the  unique  features  of  the  problem  of 
test  generation.  First  we  introduced  a  novel  on-line  method 
of  computing  the  coverage  of  the  state  space  by  the  random 
tree.  The  method  closely  approximates  known  coverage 
measures  and  is  efficient  to  compute.  Because  the  algorithm 
is  a  semi-decision  process  the  coverage  measure  is  useful  as 
termination  criteria.  In  addition  we  proposed  a  new  measure 
of  tree  growth  which  is  used  to  alter  the  search  process. 
Traditionally  RRT  methods  generate  samples  according  to 
a  uniform  distribution.  We  use  the  growth  measures  to 
adjust  the  distribution  adaptively.  As  the  growth  rate  of  the 
tree  declines  we  adjust  the  sampling  distribution  to  be  less 
biased.  This  adaptive  search  strategy  varies  bias  between 
“greedy”  and  global,  often  finding  test  trajectories  more 
quickly  than  the  traditional  algorithm. 
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